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Abstract 

The effects of free surfaces on the dynamic critical behavior of isotropic Heisenberg 
ferromagnets are studied via phenomenological scaling theory, field-theoretic renor- 
malization group tools, and high-precision computer simulations. An appropriate 
semi-infinite extension of the stochastic model J is constructed, the boundary terms 
of the associated dynamic field theory are identified, its renormalization in d < 6 
dimensions is clarified, and the boundary conditions it satisfies are given. Scaling 
laws are derived which relate the critical indices of the dynamic and static infrared 
singularities of surface quantities to familiar static bulk and surface exponents. Ac- 
curate computer-simulation data are presented for the dynamic surface structure fac- 
tor; these are in conformity with the predicted scaling behavior and could be checked 
by appropriate scattering experiments. 

PACS numbers: 75.10.Hk, 68.35.Rh, 64.60.Ht, 05.70.Jk 
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A key ingredient of modern theories of critical phenomena is the arrangement of microscop- 
ically different systems in universality classes [Jl]]. First developed for static equilibrium bulk 
critical behavior, this classification scheme has subsequently been extended both to dynamic bulk 
critical behavior [@] as well as to static surface critical behavior of semi-infinite systems at bulk 
critical points (T5],[fp . 

Since distinct dynamics may have the same equilibrium distribution, each static bulk univer- 
sality class generally splits up into several dynamic ones, represented by stochastic models called 
A, B, . . . , J [|2p. Likewise, to which static surface universality class a given system belongs is 
decided by its static bulk universality class and additional relevant surface properties. Hence each 
static surface universality class and each dynamic bulk universality class generally splits up into 
separate dynamic surface universality classes 

Unfortunately, detailed investigations of dynamic surface critical behavior have remained 
scarce and largely focused on models with purely relaxational dynamics. Isotropic Heisen- 
berg ferromagnets form an important class of systems, characterized by the presence of nondis- 
sipative (mode-coupling) terms and a conserved order parameter, whose dynamic surface critical 
behavior has not yet been investigated. To fill this gap, we shall use two different lines of ap- 
proaches: (i) phenomenological scaling and the field-theoretic renormalization group (RG); (ii) 
computer- simulation studies of the dynamic surface structure function. 

Building on (i), we shall conclude that the critical indices characterizing the dynamic surface 
critical behavior of such systems (which are 0(3) symmetric both in the bulk and at the surface) 
can be expressed in terms of known static bulk and surface critical exponents. Recently developed 
highly efficient spin dynamics algorithms [ p~Q|JTTt ] have enabled us to corroborate these findings 
numerically. 

In the simulations we utilized a classical isotropic Heisenberg ferromagnet on a simple cu- 
bic lattice whose sites i = (ii, i 2 , i 3 ) , ii,i 2 ,h = 0, . . . , L — 1, are occupied by spins Si = 
(S? , a—1, 2, 3) of length \Si\ = 1. Free boundary conditions apply along the i 3 direction, and 
periodic ones along the others, so that the layers i 3 = and i 3 = L — 1 are free surfaces. The 
dynamics is defined through 

dSi dHi at 

with the Hamiltonian 

7Yi at = — J Si-Sj — Ji Si-Sj , (2) 

<i,J> 

whose nearest-neighbor bulk and surface bonds J and J x (measured in units of temperature T) are 
ferromagnetic. 

In the thermodynamic limit the model undergoes a continuous bulk phase transition whose 
critical behavior is representative of the 0(3) universality class. Owing to the 0(3) symmetry, the 
surface of such a d=3-dimensional system cannot spontaneously order for J\j J < oo. Hence the 
surface transition that occurs at the bulk critical point T = T c is the so-called ordinary one [f|]. Its 
critical indices can be expressed in terms of two independent bulk exponents, e.g., r] and u, and 
one surface exponent, e.g., the correlation exponent ^° rd . 

In our computer simulations, ([I]) is integrated numerically for a given set of more than 700 
initial spin configurations generated by a Monte-Carlo simulation of the thermal equilibrium state 
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<7* (r; z, z'; \t-t'\) = (S»(f) Sj(f))° nn (3) 

with r = (ix— z 2 , &x — ^2)' 2 = *3> an d -z' — ^3 is calculated, where t and t' are times to which the 
initial spin configuration at t — has evolved according to ([]]). The average (.) is taken over the 
set of initial configurations. Before turning to the results, let us discuss scaling. 

Consider a d-dimensional analog of the above model in the thermodynamic limit L — > 00. Ac- 
cording to scaling considerations [ |Tl| , [Flx] and RG work [ jI3|JI5p , the dynamic bulk critical exponent 
3 is given by 

(d + 2-rj)/2 if d < d* = 6. (4) 
4 if d> 6. 

Here dj = 6 is the upper critical dimension of model J, and the mean-field (MF) value 77 = 
applies for d > d* t = 4, the static upper critical dimension. 

Let 2 < d < 4, so that hyperscaling holds. A scaling ansatz for the cumulant (Q) at T c must 
comply with the scaling forms of (a) the bulk cumulant that @ approaches as z, z' — > 00 with z—z' 
fixed, and (b) the static cumulant C(r; z, z'\ 0). Both requirements impose strong restrictions. To 
see this, recall that the static critical behavior is captured by the semi-infinite 3-vector model with 
Hamiltonian 
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defined on the half-space W\_ = {(seii, z) G R | z > 0} with boundary plane £> at 2; = 0. In this 
model there is ordinary surface critical behavior at bulk criticality if c > c sp . At the corresponding 
ordinary fixed point c, the renormalized analog of cq— c sp , takes the value c = 00. This implies 
that the order-parameter field cj) asymptotically satisfies the Dirichlet boundary condition <f>\ B = 0. 
The behavior of <f)(x) clS X = || , Z ) approaches the surface point x B = (xn, 0) follows from the 
boundary operator expansion (BOE) [Q] 

0(:r) = J D(z)a n 0( : r B ) + ... • (6) 

The scaling dimensions of the operators <fr an d d n <p are (d—2+r])/2 and (d— 2+?7|| rd )/2, respec- 

tively, whence D(^) = Z^o -2 11 • 

From requirement (b) we draw the important conclusion that and the surface operator 
d n 4> ought to retain their scaling dimensions in the dynamic case. The reason is that for the locally 
scale invariant theory considered here the scaling dimension of any local scaling operator 0(x, t) 
cannot differ from its bulk value except for points x on the boundary [|J| or at an initial condition 
HTTp. However, the initial configurations 4>t=o are drawn from the equilibrium distribution and 
should therefore not give rise to different scaling dimensions. 

Requirement (a) tells us that the dependence on t must involve the usual bulk scaling variables, 
such as tr~ s . Additional relevant or marginal dynamic surface scaling fields can be ruled out 
because (i) at the ordinary transition the surface does not act as an independent source of critical 
behavior [^] and (ii) the 0(3) invariance of the interactions precludes a local violation of the 
order-parameter conservation at the surface. Hence, sufficiently close to T c , possible times scales 
generated by the surface should be much smaller than the characteristic bulk time scale. An 
appropriate scaling ansatz for C a>3 = 5 a pC at T = T c thus is 
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C(r; z, z'\ t) = r 2 - d ^ T (z/r, z'/r; tr~ 3 ) . (7) 

The BOE (g) leads to T (z, z'; t) w (zz') ( ' ? r~' 7)/2 T (t) in the limit z, z' -> 0. By consistency we 
must have 

d n (p,u)=p<- 1 - i <T(up->) , (8) 

Cu (0, u) = const ( J+1_1, f) /3 , (9) 

where Cn(p, cu), the dynamic surface structure function, is the Fourier transform of C(r; 0, 0; t) 
with respect to r and t. 

For 4 < d < 6, the breakdown of hyperscaling (caused by the dangerous character of the 
irrelevant scaling field u ~ u ) must be taken into account. While static exponents like 77 and 77° rd 
take their MF values and 2, respectively, the value (d+2)/2 of 3 agrees with the MF value 3 of 
the magnetic shift exponent only for d = 4. 

To put the above findings on a firmer basis, let us construct an appropriate semi-infinite exten- 
sion of model J. We may assume that the surface-induced modifications of both the interactions 
and the dynamics are restricted to the immediate vicinity of the boundary B. Hence, for points 
with z > 0, we use the stochastic bulk equation [|T8"|] 



0(*,t)-C = A o (A^ + /o^x<^ , (10) 

where £ is a Gaussian random force with average (£) = and 

(C(x, t) C p (x', t')) = -A S a(3 A 5{x-x') 8{t-t') . (11) 

The derivative involves a contribution 5(z)(c — d z )4> implied by the boundary term of SH, 
which yields one obvious surface contribution to <JT0|) . But there may be others, corresponding to 
local changes of the dynamics. As expounded in [j^], the requirements of detailed balance, locality, 
order-parameter conservation, absence of irrelevant and redundant operators, and here also 0(3) 
symmetry, impose strong constraints, which are best dealt with on the level of the equivalent path- 
integral formulation (see, e.g., []i~6|]). To ensure detailed balance, the action must have the form 



J = dt U 



1 ^ fsn i\ sn 



(12) 



where J x comprises both volume and surface integrals, is the usual auxiliary field needed in 
such a Lagrangian formulation, and a prepoint discretization in time is used. 
The right-hand side of ( fTTJ| ) can be written as —IZ ■ with 

U a(3 = -A (5 Q/3 A + f e a ^ 7 ) . (13) 

To find out which surface terms must be included in J and the associated boundary conditions, we 
proceed as in Refs. As reaction operator in the action ([T2|) we use (JT3|), with the Laplacian 

A replaced by — V-V (where V acts to the left). Contributions to TZ al3 /A of the form c * l3 5(z) 
correspond to nonconservative dynamic surface terms [p and are ruled out by the presumed spin 
isotropy. As boundary conditions (valid in an operator sense [|],|5|]) we obtain 
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(<9 n -c o )0 = O, (14) 

= dn ( To + f 10,2 " A )^ = ° ' ° 5) 

<9„0 = O, (16) 
(d n - c )j>-n = A (c - d n ) (A0 - f 0x0) = , (17) 

where d n (=d z ) is a derivative along the inner normal. The first, (JT^|), is known from statics. 
The second, ([15]), means that nj", the normal component of the current j a = — A (Vj?£ + 
/o e a ^ 7 0^V(/» 7 ), whose negative divergence gives the right-hand side of (TO), vanishes at £>. [The 
precession term yields a contribution oc e a ^ (f) 13 d n cp which is zero by ([14])]. The remaining two 
ensure self-adjointness of A and consistency with the fluctuation-dissipation theorem (FDT) 

-9(t)(^(x,t)^(x',0)) = (r(x,t)(^lZf(x',0)). (18) 



To study the ordinary transition, it is convenient to set c =oo. Then (|14J) and ( |17| ) simplify 
to Dirichlet boundary conditions for and 0-7?.. In applying the RG, we consider two distinct 
renormalization schemes: a massless one, RSi, based on the e expansion about dj, and a massive 
one, RS 2 , in fixed dimension d < 4. Common to both is the form of the reparametrizations = 

Zl /2 (f> K n, = Z~ 1,2 4> Kn , d n <j> = [Z^^^U, A = ^ 4 Z X X, and / = /j,^ 2 Z 1 / 2 Z^ 1 f , 
where ji is either an arbitrary momentum scale (RSi) or the renormalized mass (RS2). That the 
Z-factors of and can be chosen reciprocal to each other follows from the fact that the bulk 
vertex function T*~\ has no primitive uv divergence oc 5'(t—i). The form of the Z-factor of / 
is a consequence of dl8| ) [|T6|]; it implies that the beta function /?/ = //d M | / becomes /?/ 



6— r)<f,+2r)\)f, with 77^ = ^^| lnZ^A. From the resulting exact relation r/l = |(6— d+r]^) 
among the values 77^ = rj and r/^ = 4 — 3 at the infrared- stable fixed point, expression (Q) for 3 
follows. 

Consider first RSi for 4 < d < 6. For u = 0, the theory is renormalizable and Z^—Zx^—l 
to all orders in /. Only one bulk renormalization factor Z\ remains whose expansion to order f 2 
may be gleaned from [|T6|], which also tells us that single insertions of the composite operator 0x 
appearing in 0-7?. additionally require a subtraction oc A0. Together with this latter counterterm, 
the above reparametrizations suffice to cure the uv singularities of the correlation and response 
functions [|T9|], including those involving the surface operators d n <fi and 0| 2= o and single insertions 
of 0x 0. The proof [ |20| ] utilizes power counting, the boundary conditions (|T4[)-([T7|), and the FDT 
[grj. The BOE (§) holds with D(z) ~ z. 

If n 7^ and d < 4, perturbative massless RG schemes are plagued by infrared singularities 



Q22[ ] . The usual escape via an e expansion here is not possible because of the different upper critical 
dimensions 6 and 4 associated with f and uq. We therefore use RS2. Writing m = /i d_4 Z u M, we 
fix all static bulk and surface counterterms (mass shift r — /i 2 , Z^, Z u , Zi t00 ) as in [ |23| , Eqs. (3.3a- 
d), (710a,b)]. The remaining dynamic (bulk) counterterms (Z\, subtraction for 0x0) can be fixed 
by appropriate (massive) normalization conditions for the dynamic bulk theory. From the resulting 
RG equations of the dynamic theory and well-established RG results for the static theory the BOE 
@ and the scaling forms (Q)-® follow. 

Our simulation results are displayed in Figs. [I] and^. They were obtained for T = T c , L = 60, 
and a total integration time of about 5000 / J. The algorithm was parallelized and implemented 
on the ALiCE cluster at the BUGH Wuppertal, where the MPI (message passing interface) library 
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was used for communication between the processes. Finite-size effects turned out to be negligible 
for this choice of parameters. 
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FIG. 1. Structure function Ch(0,cj) for J x /J = 0.3 (x), 0.73 (+), and 1.0 (*). Error bars (one 
standard deviation) are smaller than the symbol sizes. The solid lines indicate the theoretically expected 
power law (ph for uj — > 0. 



In Fig. |l] the structure function Cn(0,u}) is shown for different values of Ji/J in com- 
parison with @. For the exponent (3+I— ^|| rd )/3, we used the value 0.856 ±0.005 that fol- 
lows from the estimate $(d=3) = 2.482 ±0.002, obtained by substitution of the value p4] ] 
rj(d=3) = 0.036 ± 0.004 into (@), and the current estimate [Q r]^ ld (d=3) = 1.358 ± 0.012. 

Depending on whether Ji/J is small (Ji/J=0.3, x) or larger (Ji/J=l, *), the approach to 
the asymptotic power law @ is from above or below. In the latter case, the asymptotic regime 
is not reached within the frequency range shown in Fig. []]. The best agreement with (Q) over the 
largest frequency range is obtained for the intermediate value J\j J = 0.73 (+). 

Fig. ^ shows a scaling plot of C n (p, w) for p = (||, 0, 0). The scaling regime in x shrinks as 
the mode index n is increased from l(x) to 4(D). The shape of the scaling function in (§) , for 
x < 1, is captured surprisingly well by 



a(x) = (T [l ± (x/x ) 4 ] 



(19) 



whose form is inspired by the known zero-loop result IQgp. The exponent at the square bracket is 
chosen so as to reproduce @ in the limit x — ► 00 (p — > at fixed u 7^ 0). The amplitude a and 
the crossover parameter x are used as fit parameters. The fit shown in Fig. || has been obtained 
from the data for n = 3 and x < 1. 
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FIG. 2. Scaling function a(x) according to Data obtained for J\/J = 0.73 and p = (^,0,0), 
with n = 1, . . . , 4, are shown. Error bars (one standard deviation) are smaller than the symbol sizes. The 
solid line displays a fit to (|l^). The data for x > 1 are outside the scaling regime. 



The agreement between the data displayed in Figs. |T] and ^| and the scaling laws (§) and (§) 
is quite satisfactory, albeit small deviations are seen to remain on closer inspection. Just as in 
[|T2|], intermediate values of J\j J (= 0.73) appear to be empirically optimal in that they yield the 
largest asymptotic scaling regime. Possible sources of deviations are (i) insufficient momentum 
resolution, (ii) insufficient frequency resolution, and (iii) corrections to scaling. Momentum reso- 
lution is intimately linked to the system size L, which despite formidable progress in simulation 
techniques, still is a serious limiting factor. Frequency resolution is limited by the total integration 
time. From inspection of Cn(p, t) (not shown) we conclude that the integration time is sufficiently 
long. The frequency resolution Sco/J ~ 1.2x 1CT 3 available here rivals that of neutron scattering 
experiments JlTj]. Momentum resolution is more modest. 

Corrections to scaling may have familiar static roots or be of genuine dynamic origin. Exam- 
ples of the former category are corrections ~ (u— u*) governed by the static Wegner exponent 
uj u = f3' u (u*) ~ 0.8, for d = 3 [ffij], and corrections due to the finiteness of c ('nonzero extrap- 
olation length') [0]. An example of the second category are Wegner-type corrections ~ (/—/*), 
governed by Uf = {3f(u*, /*), the dynamic analog of u u , whose e expansion about d—6 reads 
Uf = e + 0(e 2 ) [jTfJ. While we cannot rule out that corrections of the latter category are large, we 
are not aware of compelling reasons to expect this. Thus the mentioned corrections of static origin 
may well be the most important ones. 
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In summary, we have presented a detailed study of the surface critical behavior of isotropic 
Heisenberg ferromagnets based on phenomenological scaling and field-theoretical RG methods 
and corroborated our findings through high-precision simulations. Our results depicted in Figs. [1] 
and ^ may serve as guidelines for careful experimental tests. For these, one must choose ferro- 
magnets for which dipolar forces (ignored here) can be trusted to be unimportant. 

We thank K. Wiese for helpful discussions. M.K. is indebted to the Institut fur angewandte 
Informatik at the BUGH Wuppertal for providing access to the parallel cluster ALiCE. Partial 
support by DFG (for M.K. via the Heisenberg program under grant # Kr 1322/2-1, for H.K. and 
H.W.D. via the Leibniz program and SFB 237) is gratefully acknowledged. 
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